home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 / Ham Radio 2000.iso / ham2000 / misc / dspice0s / dcdcmp.c < prev    next >
C/C++ Source or Header  |  1992-11-21  |  24KB  |  793 lines
  1. /* dcdcmp.f -- translated by f2c (version of 3 February 1990  3:36:42).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "f2c.h"
  7.  
  8. /* Common Block Declarations */
  9.  
  10. struct {
  11.     integer ielmnt, isbckt, nsbckt, iunsat, nunsat, itemps, numtem, isens, 
  12.         nsens, ifour, nfour, ifield, icode, idelim, icolum, insize, 
  13.         junode, lsbkpt, numbkp, iorder, jmnode, iur, iuc, ilc, ilr, 
  14.         numoff, isr, nmoffc, iseq, iseq1, neqn, nodevs, ndiag, iswap, 
  15.         iequa, macins, lvnim1, lx0, lvn, lynl, lyu, lyl, lx1, lx2, lx3, 
  16.         lx4, lx5, lx6, lx7, ld0, ld1, ltd, imynl, imvn, lcvn, nsnod, 
  17.         nsmat, nsval, icnod, icmat, icval, loutpt, lpol, lzer, irswpf, 
  18.         irswpr, icswpf, icswpr, irpt, jcpt, irowno, jcolno, nttbr, nttar, 
  19.         lvntmp;
  20. } tabinf_;
  21.  
  22. #define tabinf_1 tabinf_
  23.  
  24. struct {
  25.     doublereal atime, aprog[3], adate, atitle[10], defl, defw, defad, defas, 
  26.         rstats[50];
  27.     integer iwidth, lwidth, nopage;
  28. } miscel_;
  29.  
  30. #define miscel_1 miscel_
  31.  
  32. struct {
  33.     integer locate[50], jelcnt[50], nunods, ncnods, numnod, nstop, nut, nlt, 
  34.         nxtrm, ndist, ntlin, ibr, numvs, numalt, numcyc;
  35. } cirdat_;
  36.  
  37. #define cirdat_1 cirdat_
  38.  
  39. struct {
  40.     doublereal omega, time, delta, delold[7], ag[7], vt, xni, egfet, xmu, 
  41.         sfactr;
  42.     integer mode, modedc, icalc, initf, method, iord, maxord, noncon, iterno, 
  43.         itemno, nosolv, modac, ipiv, ivmflg, ipostp, iscrch, iofile;
  44. } status_;
  45.  
  46. #define status_1 status_
  47.  
  48. struct {
  49.     integer iprnta, iprntl, iprntm, iprntn, iprnto, limtim, limpts, lvlcod, 
  50.         lvltim, itl1, itl2, itl3, itl4, itl5, itl6, igoof, nogo, keof;
  51. } flags_;
  52.  
  53. #define flags_1 flags_
  54.  
  55. struct {
  56.     doublereal twopi, xlog2, xlog10, root2, rad, boltz, charge, ctok, gmin, 
  57.         reltol, abstol, vntol, trtol, chgtol, eps0, epssil, epsox, pivtol,
  58.          pivrel;
  59. } knstnt_;
  60.  
  61. #define knstnt_1 knstnt_
  62.  
  63. struct {
  64.     integer idebug[20];
  65. } debug_;
  66.  
  67. #define debug_1 debug_
  68.  
  69. struct {
  70.     doublereal value[200000];
  71. } blank_;
  72.  
  73. #define blank_1 blank_
  74.  
  75. /* Table of constant values */
  76.  
  77. static integer c__1 = 1;
  78.  
  79. /*<       subroutine dcdcmp >*/
  80. /* Subroutine */ int dcdcmp_()
  81. {
  82.     /* Format strings */
  83.     static char fmt_51[] = "(\0020*error*:  maximum entry in this column at \
  84. step \002,i4,\002 (\002,1pd12.6,\002) is less than pivtol\002)";
  85.     static char fmt_92[] = "(\0020*abort*:  pivot not in dcdcmp\002)";
  86.     static char fmt_221[] = "(\002 pivot change on fly: n= \002,i5,\002 nxti\
  87. = \002,i5,\002 nxtj= \002,i5,\002 iterno= \002,i5,\002 time= \002,1pd12.5)";
  88.  
  89.     /* System generated locals */
  90.     integer i_1;
  91.     doublereal d_1;
  92.  
  93.     /* Builtin functions */
  94.     integer s_wsfe(), do_fio(), e_wsfe();
  95.  
  96.     /* Local variables */
  97.     static integer load, locc, loci, noff, locr, iops;
  98.     static doublereal vmax;
  99.     static integer nxti, nxtj, i, j, n, noffc, locij, ifill, locnn, noffr, 
  100.         minop, i1;
  101.     static doublereal t1;
  102.     static integer j1, ispot;
  103.     extern integer indxx_();
  104.     static integer n1, n2;
  105.     static doublereal t2;
  106.     static integer isize1, isize2, lc, nstop1, no, lr;
  107. #define nodplc ((integer *)&blank_1)
  108. #define cvalue ((complex *)&blank_1)
  109.     extern /* Subroutine */ int second_(), dmpmat_();
  110.     static doublereal epsrel;
  111.     extern /* Subroutine */ int swapij_(), sizmem_(), reserv_(), extmem_();
  112.     static integer loc;
  113.     extern /* Subroutine */ int matloc_();
  114.     static doublereal perspa;
  115.     static integer nop;
  116.  
  117.     /* Fortran I/O blocks */
  118.     static cilist io__12 = { 0, 0, 0, fmt_51, 0 };
  119.     static cilist io__20 = { 0, 0, 0, fmt_92, 0 };
  120.     static cilist io__40 = { 0, 0, 0, fmt_221, 0 };
  121.  
  122.  
  123. /*<       implicit double precision (a-h,o-z) >*/
  124.  
  125. /*    this routine swaps rows and columns in the coefficient matrix 
  126. accor-*/
  127. /*ding to the numerical pivoting and minimum fillin terms.it then 
  128. performs*/
  129. /* an in-place lu factorization of the coefficient matrix. */
  130.  
  131. /* spice version 2g.6  sccsid=tabinf 3/15/83 */
  132. /*<       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, >*/
  133. /*<      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, >*/
  134. /*<      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, >*/
  135. /*<      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, >*/
  136. /*<      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, >*/
  137. /*<      5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval, >*/
  138. /*<      6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt, >*/
  139. /*<      7   irowno,jcolno,nttbr,nttar,lvntmp >*/
  140. /* spice version 2g.6  sccsid=miscel 3/15/83 */
  141. /*<       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, >*/
  142. /*<      1  defas,rstats(50),iwidth,lwidth,nopage >*/
  143. /* spice version 2g.6  sccsid=cirdat 3/15/83 */
  144. /*<       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop, >*/
  145. /*<      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc >*/
  146. /* spice version 2g.6  sccsid=status 3/15/83 */
  147. /*<       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, >*/
  148. /*<      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon, >*/
  149. /*<      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile >*/
  150. /* spice version 2g.6  sccsid=flags 3/15/83 */
  151. /*<       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts, >*/
  152. /*<      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,itl6,igoof,nogo,keof >*/
  153. /* spice version 2g.6  sccsid=knstnt 3/15/83 */
  154. /*<       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, >*/
  155. /*<      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox, >*/
  156. /*<      2   pivtol,pivrel >*/
  157. /* spice version 2g.6  sccsid=debug 3/15/83 */
  158. /*<       common/debug/ idebug(20) >*/
  159. /* spice version 2g.6  sccsid=blank 3/15/83 */
  160. /*<       common /blank/ value(200000) >*/
  161. /*<       integer nodplc(64) >*/
  162. /*<       complex cvalue(32) >*/
  163. /*<       equivalence (value(1),nodplc(1),cvalue(1)) >*/
  164.  
  165. /*<       call second(t1) >*/
  166.     second_(&t1);
  167. /*<       if (ipiv.le.0) go to 12 >*/
  168.     if (status_1.ipiv <= 0) {
  169.     goto L12;
  170.     }
  171. /*<       if (idebug(11).le.0) go to  3 >*/
  172.     if (debug_1.idebug[10] <= 0) {
  173.     goto L3;
  174.     }
  175. /*<       call dmpmat(6hdcdcmp) >*/
  176.     dmpmat_("dcdcmp", 6L);
  177. /*<       idebug(11)=idebug(11)-1 >*/
  178.     --debug_1.idebug[10];
  179. /*<     3 do 10 i=2,nstop >*/
  180. L3:
  181.     i_1 = cirdat_1.nstop;
  182.     for (i = 2; i <= i_1; ++i) {
  183. /*<       no=0 >*/
  184.     no = 0;
  185. /*<       loc=nodplc(jcpt+i) >*/
  186.     loc = nodplc[tabinf_1.jcpt + i - 1];
  187. /*<     5 if (loc.eq.0) go to 7 >*/
  188. L5:
  189.     if (loc == 0) {
  190.         goto L7;
  191.     }
  192. /*<       no=no+1 >*/
  193.     ++no;
  194. /*<       loc=nodplc(jcpt+loc) >*/
  195.     loc = nodplc[tabinf_1.jcpt + loc - 1];
  196. /*<       go to 5 >*/
  197.     goto L5;
  198. /*<     7 nodplc(numoff+i)=no >*/
  199. L7:
  200.     nodplc[tabinf_1.numoff + i - 1] = no;
  201. /*<       no=0 >*/
  202.     no = 0;
  203. /*<       loc=nodplc(irpt+i) >*/
  204.     loc = nodplc[tabinf_1.irpt + i - 1];
  205. /*<     8 if (loc.eq.0) go to 9 >*/
  206. L8:
  207.     if (loc == 0) {
  208.         goto L9;
  209.     }
  210. /*<       no=no+1 >*/
  211.     ++no;
  212. /*<       loc=nodplc(irpt+loc) >*/
  213.     loc = nodplc[tabinf_1.irpt + loc - 1];
  214. /*<       go to 8 >*/
  215.     goto L8;
  216. /*<     9 nodplc(nmoffc+i)=no >*/
  217. L9:
  218.     nodplc[tabinf_1.nmoffc + i - 1] = no;
  219. /*<    10 continue >*/
  220. /* L10: */
  221.     }
  222. /*<    12 n=1 >*/
  223. L12:
  224.     n = 1;
  225.  
  226. /*     find next pivot */
  227.  
  228. /*<    15 n=n+1 >*/
  229. L15:
  230.     ++n;
  231. /*<       if (ipiv.gt.0) go to 20 >*/
  232.     if (status_1.ipiv > 0) {
  233.     goto L20;
  234.     }
  235. /*<       if (idebug(13).le.0) go to 17 >*/
  236.     if (debug_1.idebug[12] <= 0) {
  237.     goto L17;
  238.     }
  239. /*<       call dmpmat(6hdcdcm2) >*/
  240.     dmpmat_("dcdcm2", 6L);
  241. /*<       idebug(13)=idebug(13)-1 >*/
  242.     --debug_1.idebug[12];
  243. /*<    17 if (idebug(14).le.0) go to 18 >*/
  244. L17:
  245.     if (debug_1.idebug[13] <= 0) {
  246.     goto L18;
  247.     }
  248. /*<       if (mode.ne.2) go to 18 >*/
  249.     if (status_1.mode != 2) {
  250.     goto L18;
  251.     }
  252. /*<       call dmpmat(6hdcdcm3) >*/
  253.     dmpmat_("dcdcm3", 6L);
  254. /*<       idebug(14)=idebug(14)-1 >*/
  255.     --debug_1.idebug[13];
  256. /*<    18 if (idebug(15).le.0.or.icalc.le.10) go to 19 >*/
  257. L18:
  258.     if (debug_1.idebug[14] <= 0 || status_1.icalc <= 10) {
  259.     goto L19;
  260.     }
  261. /*<       call dmpmat(6hdcdcm4) >*/
  262.     dmpmat_("dcdcm4", 6L);
  263. /*<       idebug(15)=idebug(15)-1 >*/
  264.     --debug_1.idebug[14];
  265. /*<    19 nxti=n >*/
  266. L19:
  267.     nxti = n;
  268. /*<       nxtj=n >*/
  269.     nxtj = n;
  270. /*<       go to 120 >*/
  271.     goto L120;
  272.  
  273. /*     search the coresponding column for max entry */
  274.  
  275. /*<    20 vmax=0.0d0 >*/
  276. L20:
  277.     vmax = 0.;
  278. /*<       loci=n >*/
  279.     loci = n;
  280. /*<    25 loci=nodplc(irpt+loci) >*/
  281. L25:
  282.     loci = nodplc[tabinf_1.irpt + loci - 1];
  283. /*<       if (loci.eq.0) go to 50 >*/
  284.     if (loci == 0) {
  285.     goto L50;
  286.     }
  287. /*<       i=nodplc(irowno+loci) >*/
  288.     i = nodplc[tabinf_1.irowno + loci - 1];
  289. /*<       if (i.lt.n) go to 25 >*/
  290.     if (i < n) {
  291.     goto L25;
  292.     }
  293. /*<    30 if (dabs(value(lvn+loci)).le.vmax) go to 25 >*/
  294. /* L30: */
  295.     if ((d_1 = blank_1.value[tabinf_1.lvn + loci - 1], abs(d_1)) <= vmax) {
  296.     goto L25;
  297.     }
  298. /*<       vmax=dabs(value(lvn+loci)) >*/
  299.     vmax = (d_1 = blank_1.value[tabinf_1.lvn + loci - 1], abs(d_1));
  300. /*<       go to 25 >*/
  301.     goto L25;
  302. /*<    50 if (vmax.gt.pivtol) go to 60 >*/
  303. L50:
  304.     if (vmax > knstnt_1.pivtol) {
  305.     goto L60;
  306.     }
  307. /*<       write(iofile,51) n,vmax >*/
  308.     io__12.ciunit = status_1.iofile;
  309.     s_wsfe(&io__12);
  310.     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
  311.     do_fio(&c__1, (char *)&vmax, (ftnlen)sizeof(doublereal));
  312.     e_wsfe();
  313. /*<    51 format('0*error*:  maximum entry in this column at step ',i4,' (', >*/
  314. /*<      1   1pd12.6,') is less than pivtol') >*/
  315. /*<       igoof=1 >*/
  316.     flags_1.igoof = 1;
  317. /*<       return >*/
  318.     return 0;
  319. /*<    60 epsrel=dmax1(pivrel*vmax,pivtol) >*/
  320. L60:
  321. /* Computing MAX */
  322.     d_1 = knstnt_1.pivrel * vmax;
  323.     epsrel = max(knstnt_1.pivtol,d_1);
  324. /*<       if (n.ge.nstop) go to 200 >*/
  325.     if (n >= cirdat_1.nstop) {
  326.     goto L200;
  327.     }
  328. /*<       if (ipiv.le.0) go to 120 >*/
  329.     if (status_1.ipiv <= 0) {
  330.     goto L120;
  331.     }
  332.  
  333. /*     pivoting on the diagonal */
  334.  
  335. /*<       minop=100000 >*/
  336.     minop = 100000;
  337. /*<       nxti=0 >*/
  338.     nxti = 0;
  339. /*<       do 70 i=n,nstop >*/
  340.     i_1 = cirdat_1.nstop;
  341.     for (i = n; i <= i_1; ++i) {
  342. /*<       i1=nodplc(irswpf+i) >*/
  343.     i1 = nodplc[tabinf_1.irswpf + i - 1];
  344. /*<       j1=nodplc(icswpf+i) >*/
  345.     j1 = nodplc[tabinf_1.icswpf + i - 1];
  346. /*<       ispot=indxx(i1,j1) >*/
  347.     ispot = indxx_(&i1, &j1);
  348. /*<       if (ispot.eq.1) go to 70 >*/
  349.     if (ispot == 1) {
  350.         goto L70;
  351.     }
  352. /*<       if (dabs(value(lvn+ispot)).lt.epsrel) go to 70 >*/
  353.     if ((d_1 = blank_1.value[tabinf_1.lvn + ispot - 1], abs(d_1)) < 
  354.         epsrel) {
  355.         goto L70;
  356.     }
  357. /*<       nop=(nodplc(numoff+i)-1)*(nodplc(nmoffc+i)-1) >*/
  358.     nop = (nodplc[tabinf_1.numoff + i - 1] - 1) * (nodplc[tabinf_1.nmoffc 
  359.         + i - 1] - 1);
  360. /*<       if (nop.ge.minop) go to 70 >*/
  361.     if (nop >= minop) {
  362.         goto L70;
  363.     }
  364. /*<       minop=nop >*/
  365.     minop = nop;
  366. /*<       nxti=i >*/
  367.     nxti = i;
  368. /*<       nxtj=i >*/
  369.     nxtj = i;
  370. /*<       if (minop.le.0) go to 95 >*/
  371.     if (minop <= 0) {
  372.         goto L95;
  373.     }
  374. /*<    70 continue >*/
  375. L70:
  376.     ;}
  377. /*<       if (nxti.le.0) go to 75 >*/
  378.     if (nxti <= 0) {
  379.     goto L75;
  380.     }
  381. /*<       if (nxti-n) 120,120,100 >*/
  382.     if (nxti - n <= 0) {
  383.     goto L120;
  384.     } else {
  385.     goto L100;
  386.     }
  387.  
  388. /*     pivoting on the entire matrix */
  389.  
  390. /*<    75 do 90 i=n,nstop >*/
  391. L75:
  392.     i_1 = cirdat_1.nstop;
  393.     for (i = n; i <= i_1; ++i) {
  394. /*<       loc=i >*/
  395.     loc = i;
  396. /*<    80 loc=nodplc(jcpt+loc) >*/
  397. L80:
  398.     loc = nodplc[tabinf_1.jcpt + loc - 1];
  399. /*<       if (loc.eq.0) go to 90 >*/
  400.     if (loc == 0) {
  401.         goto L90;
  402.     }
  403. /*<       j=nodplc(jcolno+loc) >*/
  404.     j = nodplc[tabinf_1.jcolno + loc - 1];
  405. /*<       if (j.lt.n) go to 80 >*/
  406.     if (j < n) {
  407.         goto L80;
  408.     }
  409. /*<       if (dabs(value(lvn+loc)).lt.epsrel) go to 80 >*/
  410.     if ((d_1 = blank_1.value[tabinf_1.lvn + loc - 1], abs(d_1)) < epsrel) 
  411.         {
  412.         goto L80;
  413.     }
  414. /*<       nop=(nodplc(numoff+i)-1)*(nodplc(nmoffc+j)-1) >*/
  415.     nop = (nodplc[tabinf_1.numoff + i - 1] - 1) * (nodplc[tabinf_1.nmoffc 
  416.         + j - 1] - 1);
  417. /*<       if (nop.ge.minop) go to 80 >*/
  418.     if (nop >= minop) {
  419.         goto L80;
  420.     }
  421. /*<       minop=nop >*/
  422.     minop = nop;
  423. /*<       nxti=i >*/
  424.     nxti = i;
  425. /*<       nxtj=j >*/
  426.     nxtj = j;
  427. /*<       if (minop.le.0) go to 95 >*/
  428.     if (minop <= 0) {
  429.         goto L95;
  430.     }
  431. /*<    90 continue >*/
  432. L90:
  433.     ;}
  434. /*<       if (nxti.gt.0) go to 95 >*/
  435.     if (nxti > 0) {
  436.     goto L95;
  437.     }
  438. /*<       write (iofile,92) >*/
  439.     io__20.ciunit = status_1.iofile;
  440.     s_wsfe(&io__20);
  441.     e_wsfe();
  442. /*<    92 format('0*abort*:  pivot not in dcdcmp') >*/
  443. /*<       igoof=1 >*/
  444.     flags_1.igoof = 1;
  445. /*<       go to 200 >*/
  446.     goto L200;
  447. /*<    95 if (nxti.eq.n.and.nxtj.eq.n) go to 120 >*/
  448. L95:
  449.     if (nxti == n && nxtj == n) {
  450.     goto L120;
  451.     }
  452. /*<       if (nxti.eq.n) go to 105 >*/
  453.     if (nxti == n) {
  454.     goto L105;
  455.     }
  456.  
  457. /*     a pivot has been found */
  458.  
  459. /*<   100 load=nodplc(irswpf+nxti) >*/
  460. L100:
  461.     load = nodplc[tabinf_1.irswpf + nxti - 1];
  462. /*<       lr=nodplc(irswpf+n) >*/
  463.     lr = nodplc[tabinf_1.irswpf + n - 1];
  464. /*<       nodplc(irswpf+nxti)=lr >*/
  465.     nodplc[tabinf_1.irswpf + nxti - 1] = lr;
  466. /*<       nodplc(irswpr+lr)=nxti >*/
  467.     nodplc[tabinf_1.irswpr + lr - 1] = nxti;
  468. /*<       nodplc(irswpf+n)=load >*/
  469.     nodplc[tabinf_1.irswpf + n - 1] = load;
  470. /*<       nodplc(irswpr+load)=n >*/
  471.     nodplc[tabinf_1.irswpr + load - 1] = n;
  472. /*<       noff=nodplc(numoff+nxti) >*/
  473.     noff = nodplc[tabinf_1.numoff + nxti - 1];
  474. /*<       nodplc(numoff+nxti)=nodplc(numoff+n) >*/
  475.     nodplc[tabinf_1.numoff + nxti - 1] = nodplc[tabinf_1.numoff + n - 1];
  476. /*<       nodplc(numoff+n)=noff >*/
  477.     nodplc[tabinf_1.numoff + n - 1] = noff;
  478. /*<       if (nxtj.eq.n) go to 110 >*/
  479.     if (nxtj == n) {
  480.     goto L110;
  481.     }
  482. /*<   105 load=nodplc(icswpf+nxtj) >*/
  483. L105:
  484.     load = nodplc[tabinf_1.icswpf + nxtj - 1];
  485. /*<       lc=nodplc(icswpf+n) >*/
  486.     lc = nodplc[tabinf_1.icswpf + n - 1];
  487. /*<       nodplc(icswpf+nxtj)=lc >*/
  488.     nodplc[tabinf_1.icswpf + nxtj - 1] = lc;
  489. /*<       nodplc(icswpr+lc)=nxtj >*/
  490.     nodplc[tabinf_1.icswpr + lc - 1] = nxtj;
  491. /*<       nodplc(icswpf+n)=load >*/
  492.     nodplc[tabinf_1.icswpf + n - 1] = load;
  493. /*<       nodplc(icswpr+load)=n >*/
  494.     nodplc[tabinf_1.icswpr + load - 1] = n;
  495. /*<       noff=nodplc(nmoffc+nxtj) >*/
  496.     noff = nodplc[tabinf_1.nmoffc + nxtj - 1];
  497. /*<       nodplc(nmoffc+nxtj)=nodplc(nmoffc+n) >*/
  498.     nodplc[tabinf_1.nmoffc + nxtj - 1] = nodplc[tabinf_1.nmoffc + n - 1];
  499. /*<       nodplc(nmoffc+n)=noff >*/
  500.     nodplc[tabinf_1.nmoffc + n - 1] = noff;
  501. /*<   110 call swapij(nxti,n,nxtj,n) >*/
  502. L110:
  503.     swapij_(&nxti, &n, &nxtj, &n);
  504. /*<       nxti=n >*/
  505.     nxti = n;
  506. /*<       nxtj=n >*/
  507.     nxtj = n;
  508.  
  509. /*     calculate contribution from nxti, nxtj and find fill-ins */
  510.  
  511. /*<   120 if (n.ge.nstop) go to 200 >*/
  512. L120:
  513.     if (n >= cirdat_1.nstop) {
  514.     goto L200;
  515.     }
  516. /*<       n1=nodplc(irswpf+nxti) >*/
  517.     n1 = nodplc[tabinf_1.irswpf + nxti - 1];
  518. /*<       n2=nodplc(icswpf+nxtj) >*/
  519.     n2 = nodplc[tabinf_1.icswpf + nxtj - 1];
  520. /*<       locnn=indxx(n1,n2) >*/
  521.     locnn = indxx_(&n1, &n2);
  522. /*<       if (ipiv.le.0 .and. dabs(value(lvn+locnn)).lt.pivtol) go to 220 >*/
  523.     if (status_1.ipiv <= 0 && (d_1 = blank_1.value[tabinf_1.lvn + locnn - 1], 
  524.         abs(d_1)) < knstnt_1.pivtol) {
  525.     goto L220;
  526.     }
  527.  
  528. /*     down col j */
  529.  
  530. /*<       locr=nodplc(irpt+locnn) >*/
  531.     locr = nodplc[tabinf_1.irpt + locnn - 1];
  532. /*<   125 if (locr.eq.0) go to 180 >*/
  533. L125:
  534.     if (locr == 0) {
  535.     goto L180;
  536.     }
  537. /*<       i=nodplc(irowno+locr) >*/
  538.     i = nodplc[tabinf_1.irowno + locr - 1];
  539. /*<       value(lvn+locr)=value(lvn+locr)/value(lvn+locnn) >*/
  540.     blank_1.value[tabinf_1.lvn + locr - 1] /= blank_1.value[tabinf_1.lvn + 
  541.         locnn - 1];
  542. /*<       locc=nodplc(jcpt+locnn) >*/
  543.     locc = nodplc[tabinf_1.jcpt + locnn - 1];
  544.  
  545. /*     for each column element look up row nxti */
  546.  
  547. /*<   130 if (locc.eq.0) go to 170 >*/
  548. L130:
  549.     if (locc == 0) {
  550.     goto L170;
  551.     }
  552. /*<       j=nodplc(jcolno+locc) >*/
  553.     j = nodplc[tabinf_1.jcolno + locc - 1];
  554.  
  555. /*     check for fill-in (i,j) */
  556.  
  557. /*<       if (ipiv.le.0) go to 135 >*/
  558.     if (status_1.ipiv <= 0) {
  559.     goto L135;
  560.     }
  561. /*<       call sizmem(jcpt,isize1) >*/
  562.     sizmem_(&tabinf_1.jcpt, &isize1);
  563. /*<       call reserv(i,j) >*/
  564.     reserv_(&i, &j);
  565. /*<       call sizmem(jcpt,isize2) >*/
  566.     sizmem_(&tabinf_1.jcpt, &isize2);
  567. /*<       if (isize1.eq.isize2) go to 135 >*/
  568.     if (isize1 == isize2) {
  569.     goto L135;
  570.     }
  571. /*<       call extmem(lvn,1) >*/
  572.     extmem_(&tabinf_1.lvn, &c__1);
  573. /*<       nttbr=nttbr+1 >*/
  574.     ++tabinf_1.nttbr;
  575. /*<       value(lvn+nstop+nttbr)=0.0d0 >*/
  576.     blank_1.value[tabinf_1.lvn + cirdat_1.nstop + tabinf_1.nttbr - 1] = 0.;
  577.  
  578. /*     locate element (i,j) */
  579.  
  580. /*<   135 if (j.lt.i) go to 145 >*/
  581. L135:
  582.     if (j < i) {
  583.     goto L145;
  584.     }
  585. /*<       locij=locc >*/
  586.     locij = locc;
  587. /*<   140 locij=nodplc(irpt+locij) >*/
  588. L140:
  589.     locij = nodplc[tabinf_1.irpt + locij - 1];
  590. /*<       if (nodplc(irowno+locij).eq.i) go to 155 >*/
  591.     if (nodplc[tabinf_1.irowno + locij - 1] == i) {
  592.     goto L155;
  593.     }
  594. /*<       go to 140 >*/
  595.     goto L140;
  596. /*<   145 locij=locr >*/
  597. L145:
  598.     locij = locr;
  599. /*<   150 locij=nodplc(jcpt+locij) >*/
  600. L150:
  601.     locij = nodplc[tabinf_1.jcpt + locij - 1];
  602. /*<       if (nodplc(jcolno+locij).eq.j) go to 155 >*/
  603.     if (nodplc[tabinf_1.jcolno + locij - 1] == j) {
  604.     goto L155;
  605.     }
  606. /*<       go to 150 >*/
  607.     goto L150;
  608. /*<   155 value(lvn+locij)=value(lvn+locij)- >*/
  609. /*<      1                  value(lvn+locc)*value(lvn+locr) >*/
  610. L155:
  611.     blank_1.value[tabinf_1.lvn + locij - 1] -= blank_1.value[tabinf_1.lvn + 
  612.         locc - 1] * blank_1.value[tabinf_1.lvn + locr - 1];
  613. /*<   160 locc=nodplc(jcpt+locc) >*/
  614. /* L160: */
  615.     locc = nodplc[tabinf_1.jcpt + locc - 1];
  616. /*<       go to 130 >*/
  617.     goto L130;
  618. /*<   170 locr=nodplc(irpt+locr) >*/
  619. L170:
  620.     locr = nodplc[tabinf_1.irpt + locr - 1];
  621. /*<       if (ipiv.le.0) go to 125 >*/
  622.     if (status_1.ipiv <= 0) {
  623.     goto L125;
  624.     }
  625. /*<       nodplc(numoff+i)=nodplc(numoff+i)-1 >*/
  626.     --nodplc[tabinf_1.numoff + i - 1];
  627. /*<       go to 125 >*/
  628.     goto L125;
  629.  
  630. /*     reduce nmoffc for each element in col nxti */
  631.  
  632. /*<   180 if (ipiv.le.0) go to 15 >*/
  633. L180:
  634.     if (status_1.ipiv <= 0) {
  635.     goto L15;
  636.     }
  637. /*<       locc=nodplc(jcpt+locnn) >*/
  638.     locc = nodplc[tabinf_1.jcpt + locnn - 1];
  639. /*<   185 if (locc.eq.0) go to 15 >*/
  640. L185:
  641.     if (locc == 0) {
  642.     goto L15;
  643.     }
  644. /*<       j=nodplc(jcolno+locc) >*/
  645.     j = nodplc[tabinf_1.jcolno + locc - 1];
  646. /*<       nodplc(nmoffc+j)=nodplc(nmoffc+j)-1 >*/
  647.     --nodplc[tabinf_1.nmoffc + j - 1];
  648. /*<       locc=nodplc(jcpt+locc) >*/
  649.     locc = nodplc[tabinf_1.jcpt + locc - 1];
  650. /*<       go to 185 >*/
  651.     goto L185;
  652.  
  653. /*     done */
  654.  
  655. /*<   200 if (ipiv.eq.0) go to 210 >*/
  656. L200:
  657.     if (status_1.ipiv == 0) {
  658.     goto L210;
  659.     }
  660. /*<       if (idebug(17).le.0) go to 202 >*/
  661.     if (debug_1.idebug[16] <= 0) {
  662.     goto L202;
  663.     }
  664. /*<       call dmpmat(6hdcdcm5) >*/
  665.     dmpmat_("dcdcm5", 6L);
  666. /*<       idebug(17)=idebug(17)-1 >*/
  667.     --debug_1.idebug[16];
  668. /*<   202 call matloc >*/
  669. L202:
  670.     matloc_();
  671. /*<       rstats(49)=rstats(49)+1.0d0 >*/
  672.     miscel_1.rstats[48] += 1.;
  673. /*<       ipiv=0 >*/
  674.     status_1.ipiv = 0;
  675. /*<       if (lvlcod.eq.2) lvlcod=3 >*/
  676.     if (flags_1.lvlcod == 2) {
  677.     flags_1.lvlcod = 3;
  678.     }
  679. /*<       ifill=nttbr-nttar >*/
  680.     ifill = tabinf_1.nttbr - tabinf_1.nttar;
  681. /*<       perspa=100.0d0*(1.0d0-dble(nttbr)/dble(nstop*nstop)) >*/
  682.     perspa = (1. - (doublereal) tabinf_1.nttbr / (doublereal) (cirdat_1.nstop 
  683.         * cirdat_1.nstop)) * 100.;
  684.  
  685. /*  calculation of operation count (operation := `*' or `/'): */
  686.  
  687. /*     noffr := off-diagonal elements in row, not including diagonal, */
  688. /*                counting only those elements in the remainder matrix */
  689. /*     noffc := off-diagonal elements in column, not including diagonal, 
  690. */
  691. /*                counting only those elements in the remainder matrix */
  692.  
  693. /*     then we have */
  694.  
  695. /*       lu decomposition     requires sigma(2,nstop-1) {noffc + noffc*
  696. noffr}*/
  697. /*       forward substitution          sigma(2,nstop-1) {noffc + 1}   +   
  698. 1*/
  699. /*        backward substitution         sigma(2,nstop-1) {noffr} */
  700.  
  701. /*     which sums to */
  702.  
  703. /*              sigma(2,nstop-1) {noffc + noffc*noffr + (noffc+1) + noffr}
  704.  + 1*/
  705. /*         or */
  706. /*               sigma(2,nstop-1) {noffc*(noffr+2) + noffr + 1}   +   1 */
  707.  
  708.  
  709. /*<       iops=1 >*/
  710.     iops = 1;
  711. /*<       nstop1=nstop-1 >*/
  712.     nstop1 = cirdat_1.nstop - 1;
  713. /*<       do 205 i=2,nstop1 >*/
  714.     i_1 = nstop1;
  715.     for (i = 2; i <= i_1; ++i) {
  716. /*<       noffr=nodplc(numoff+i)-1 >*/
  717.     noffr = nodplc[tabinf_1.numoff + i - 1] - 1;
  718. /*<       noffc=nodplc(nmoffc+i)-1 >*/
  719.     noffc = nodplc[tabinf_1.nmoffc + i - 1] - 1;
  720. /*<       iops=iops+noffr+noffc*(noffr+2)+1 >*/
  721.     iops = iops + noffr + noffc * (noffr + 2) + 1;
  722. /*<   205 continue >*/
  723. /* L205: */
  724.     }
  725. /*<       rstats(20)=nstop >*/
  726.     miscel_1.rstats[19] = (doublereal) cirdat_1.nstop;
  727. /*<       rstats(21)=nttar >*/
  728.     miscel_1.rstats[20] = (doublereal) tabinf_1.nttar;
  729. /*<       rstats(22)=nttbr >*/
  730.     miscel_1.rstats[21] = (doublereal) tabinf_1.nttbr;
  731. /*<       rstats(23)=ifill >*/
  732.     miscel_1.rstats[22] = (doublereal) ifill;
  733. /*<       rstats(24)=0.0d0 >*/
  734.     miscel_1.rstats[23] = 0.;
  735. /*<       rstats(25)=nttbr >*/
  736.     miscel_1.rstats[24] = (doublereal) tabinf_1.nttbr;
  737. /*<       rstats(26)=iops >*/
  738.     miscel_1.rstats[25] = (doublereal) iops;
  739. /*<       rstats(27)=perspa >*/
  740.     miscel_1.rstats[26] = perspa;
  741. /*<       go to 215 >*/
  742.     goto L215;
  743. /*<   210 if (idebug(18).le.0) go to 212 >*/
  744. L210:
  745.     if (debug_1.idebug[17] <= 0) {
  746.     goto L212;
  747.     }
  748. /*<       call dmpmat(6hdcdcm6) >*/
  749.     dmpmat_("dcdcm6", 6L);
  750. /*<       idebug(18)=idebug(18)-1 >*/
  751.     --debug_1.idebug[17];
  752. /*<   212 if (idebug(19).le.0.or.icalc.le.10) go to 215 >*/
  753. L212:
  754.     if (debug_1.idebug[18] <= 0 || status_1.icalc <= 10) {
  755.     goto L215;
  756.     }
  757. /*<       call dmpmat(6hdcdcm7) >*/
  758.     dmpmat_("dcdcm7", 6L);
  759. /*<       idebug(19)=idebug(19)-1 >*/
  760.     --debug_1.idebug[18];
  761. /*<   215 call second(t2) >*/
  762. L215:
  763.     second_(&t2);
  764. /*<       rstats(50)=rstats(50)+t2-t1 >*/
  765.     miscel_1.rstats[49] = miscel_1.rstats[49] + t2 - t1;
  766. /*<        return >*/
  767.     return 0;
  768. /*<   220 ipiv=1 >*/
  769. L220:
  770.     status_1.ipiv = 1;
  771. /*<       write(iofile,221) n,nxti,nxtj,iterno,time >*/
  772.     io__40.ciunit = status_1.iofile;
  773.     s_wsfe(&io__40);
  774.     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
  775.     do_fio(&c__1, (char *)&nxti, (ftnlen)sizeof(integer));
  776.     do_fio(&c__1, (char *)&nxtj, (ftnlen)sizeof(integer));
  777.     do_fio(&c__1, (char *)&status_1.iterno, (ftnlen)sizeof(integer));
  778.     do_fio(&c__1, (char *)&status_1.time, (ftnlen)sizeof(doublereal));
  779.     e_wsfe();
  780. /*<   221 format(' pivot change on fly: n= ',i5,' nxti= ',i5,' nxtj= ', >*/
  781. /*<      1       i5,' iterno= ',i5,' time= ',1pd12.5) >*/
  782. /*<       rstats(49)=rstats(49)+1.0d0 >*/
  783.     miscel_1.rstats[48] += 1.;
  784. /*<       go to 20 >*/
  785.     goto L20;
  786. /*<       end >*/
  787. } /* dcdcmp_ */
  788.  
  789. #undef cvalue
  790. #undef nodplc
  791.  
  792.  
  793.